notebook.community

Edit and run



In [1]:

    
import time
import os
import sys

import numpy as np
import matplotlib
matplotlib.use('nbagg')
#from matplotlib import style
#style.use('ggplot')
import matplotlib.pyplot as plt

import astropy.units as u
from astropy import stats, wcs
from astropy.io import fits
from astropy.convolution import Gaussian2DKernel
from astropy.coordinates import SkyCoord

from mmtwfs.wfs import *
from mmtwfs.zernike import ZernikeVector
from mmtwfs.telescope import MMT
import poppy
import photutils



In [2]:

    
%load_ext autoreload
%autoreload 2



In [18]:

    
sys.path.append("/Users/tim/src/cwfs/python")
%cd /Users/tim/MMT/61_inch/20181219/









    



/Users/tim/MMT/61_inch/20181219



In [19]:

    
from lsst.cwfs.instrument import Instrument
from lsst.cwfs.algorithm import Algorithm
from lsst.cwfs.image import Image, readFile
import lsst.cwfs.plots as plots



In [49]:

    
r = 60
extra_whole = fits.open("rts2focus0012.fits")[1].data # focus = 3190
intra_whole = fits.open("rts2focus0017.fits")[1].data # focus = 1992
#intra = intra_whole[942-18:942+18,158-18:158+18] - np.median(intra_whole)
#extra = extra_whole[945-18:945+18,154-18:154+18] - np.median(extra_whole)
intra = intra_whole[780-r:780+r,503-r:503+r] - np.median(intra_whole)
extra = extra_whole[780-r:780+r,509-r:509+r] - np.median(extra_whole)



In [50]:

    
#plt.imshow(intra_whole[942-20:942+20,158-20:158+20], origin="lower")
#plt.imshow(extra_whole[945-20:945+20,154-20:154+20], origin="lower")
plt.imshow(intra, origin="lower")
plt.show()



In [6]:

    
plt.imshow(extra, origin="lower")
plt.show()



In [51]:

    
fits.writeto("intra.fits", intra, overwrite=True)
fits.writeto("extra.fits", extra, overwrite=True)

These images are binned 3x3 so the pixels are 42 um. Thus the pupil diameter is:



In [34]:

    
pupsize = 25.25 * 3 * 14
pupsize * u.um









    Out[34]:





$1060.5 \; \mathrm{\mu m}$

The mirror diameter is 1.54 meters and the focal ratio is 13.5. So the focal length is:



In [35]:

    
diameter = 1.54 * u.m
radius = diameter / 2.
focal_length = diameter * 13.5
focal_length









    Out[35]:





$20.79 \; \mathrm{m}$

Calculate the angle of beam convergence:



In [36]:

    
ang = np.arctan2(radius, focal_length)
ang









    Out[36]:





$0.037020116 \; \mathrm{rad}$



In [37]:

    
offset = pupsize / np.tan(ang) * u.um
offset









    Out[37]:





$28633.5 \; \mathrm{\mu m}$



In [12]:

    
def kuiper_focus(pup_rad, focus, binning=3, direction='out'):
    pupsize = 14 * binning * pup_rad * u.um
    diameter = 1.54 * u.m
    radius = diameter / 2
    focal_length = diameter * 13.5
    beam_ang = np.arctan2(radius, focal_length)
    offset = pupsize / np.tan(beam_ang) * u.um
    counts_per_um = 600. / 28633.5  # empirically determined
    correction = offset.value * counts_per_um
    if 'out' in direction:
        correction *= -1
    return focus + correction



In [13]:

    
kuiper_focus(25.25, 3190, direction='out')









    Out[13]:





2590.0

So about 14 mm of focus offset at focal plane for 300 counts of focus change which gives 20 counts/mm. For motion at the secondary, scale by the ratio of the squares of the focal ratios:



In [14]:

    
counts_per_mm = 600. / offset.to(u.mm).value
counts_per_mm_m2 = counts_per_mm / (13.5**2 / 4.**2)
counts_per_mm, counts_per_mm_m2









    Out[14]:





(20.954476400020955, 1.8396248142679577)



In [15]:

    
obscuration = 0.4096 / 1.54
obscuration









    Out[15]:





0.265974025974026



In [16]:

    
nmperrad = radius.to(u.nm).value
nmperasec = nmperrad / 206265.
nmperasec









    Out[16]:





3733.061837926938



In [52]:

    
fieldXY = [0., 0.]
I1 = Image(readFile("intra.fits"), fieldXY, Image.INTRA)
I2 = Image(readFile("extra.fits"), fieldXY, Image.EXTRA)



In [53]:

    
fig, ax = plt.subplots()
im = ax.imshow(I2.image, origin='lower', cmap='RdBu')
cbar = fig.colorbar(im)
fig.show()



In [54]:

    
kuiper = Instrument('61inch', I1.sizeinPix)



In [55]:

    
#kuiper.offset = offset.to(u.m).value  # needs to be focus offset between intra and extra positions
kuiper.offset









    Out[55]:





0.028634



In [56]:

    
algo = Algorithm('exp', kuiper, 3)



In [57]:

    
algo.runIt(kuiper, I1, I2, 'onAxis')









    



resetting images: I1 and I2
imageCoCenter: (x1,y1)=(   60.43,   60.70)

imageCoCenter: (x1,y1)=(   59.23,   59.72)

itr = 0, z4-z22
[-402.  -14.   -1.   51.  -20.   56.    4.   73.    2.    2.  -32.    4.
    2.  -10.   -2.   -3.    7.   -4.    8.]
itr = 1, z4-z22
[-468.  -10.   31.   40.  -16.   48.   11.   75.    1.    2.  -26.    2.
    3.   -8.   -4.   -3.    3.   -3.   -2.]
itr = 2, z4-z22
[-515.   -0.   43.   36.  -10.   43.   18.   74.   -0.    1.  -29.    3.
    2.   -7.   -4.   -3.    2.   -5.   -8.]
itr = 3, z4-z22
[-545.    8.   56.   33.   -6.   41.   24.   74.   -0.    1.  -28.    3.
    1.   -6.   -5.   -2.    2.   -5.  -10.]
itr = 4, z4-z22
[-578.   -2.   73.   74.   -4.   38.   19.   73.    1.    2.  -26.    5.
   -0.  -24.   -4.   -3.    1.   -6.  -12.]
itr = 5, z4-z22
[-505.    1.   66.   31.   -7.   53.   24.   78.    2.    1.  -27.    4.
    0.   -6.   -4.   -2.    2.   -4.   -8.]
itr = 6, z4-z22
[-462.  -19.   54.   54.  -10.   68.   19.   80.    3.    3.  -26.    4.
    1.  -17.   -3.   -2.    1.   -4.   -5.]
itr = 7, z4-z22
[-439.   -3.   48.   30.  -15.   67.   17.   80.    3.    1.  -26.    2.
    2.   -4.   -2.   -0.    2.   -4.   -1.]
itr = 8, z4-z22
[-423.  -17.   31.   54.  -14.   74.   13.   81.    4.    2.  -26.    5.
    2.  -16.   -2.   -0.    3.   -4.   -0.]
itr = 9, z4-z22
[-414.  -12.   29.   49.  -19.   75.    7.   80.    4.    2.  -33.    6.
    3.  -12.   -2.   -0.    4.   -6.    1.]
itr = 10, z4-z22
[-405.   -8.   25.   43.  -22.   75.    3.   80.    3.    1.  -36.    6.
    4.   -8.   -2.    0.    5.   -6.    2.]
itr = 11, z4-z22
[-393.  -25.   14.   68.  -23.   80.   -4.   81.    4.    3.  -40.   10.
    5.  -21.   -2.    0.    7.   -7.    2.]
itr = 12, z4-z22
[-387.  -20.   15.   61.  -26.   81.   -7.   80.    3.    2.  -42.   10.
    6.  -17.   -2.    1.    8.   -7.    4.]
itr = 13, z4-z22
[-380.  -17.   11.   54.  -29.   81.   -9.   80.    3.    2.  -43.   10.
    7.  -13.   -2.    1.    8.   -7.    4.]
itr = 14, z4-z22
[-374.  -15.    8.   48.  -32.   81.  -11.   80.    3.    2.  -43.    9.
    8.  -10.   -2.    1.    9.   -7.    5.]



In [26]:

    
print(algo.zer4UpNm)









    



[-373.82672972  -14.74702685    8.46076505   47.70437384  -31.85888618
   81.43029275  -11.09791341   80.40980238    2.78754921    1.69393212
  -43.08549841    9.4639366     8.34297213   -9.91783341   -2.05401921
    0.8499141     9.01395269   -6.93465254    5.21078698]



In [58]:

    
plots.plotZer(algo.zer4UpNm, 'nm')



In [59]:

    
zv = ZernikeVector()
zv.from_array(algo.zer4UpNm, modestart=4, normalized=True)
zv.denormalize()
zv









    Out[59]:





Fringe Coefficients
 Z04:                -647.5 nm 	 Defocus (2, 0)
 Z05:                -36.12 nm 	 Primary Astig at 45° (2, -2)
 Z06:                 20.72 nm 	 Primary Astig at 0° (2, 2)
 Z07:                 134.9 nm 	 Primary Y Coma (3, -1)
 Z08:                -90.11 nm 	 Primary X Coma (3, 1)
 Z09:                 230.3 nm 	 Y Trefoil (3, -3)
 Z10:                -31.39 nm 	 X Trefoil (3, 3)
 Z11:                 179.8 nm 	 Primary Spherical (4, 0)
 Z12:                 8.815 nm 	 Secondary Astigmatism at 0° (4, 2)
 Z13:                 5.357 nm 	 Secondary Astigmatism at 45° (4, -2)
 Z14:                -136.2 nm 	 X Tetrafoil (4, 4)
 Z15:                 29.93 nm 	 Y Tetrafoil (4, -4)
 Z16:                  28.9 nm 	 Secondary X Coma (5, 1)
 Z17:                -34.36 nm 	 Secondary Y Coma (5, -1)
 Z18:                -7.115 nm 	 Secondary X Trefoil (5, 3)
 Z19:                 2.944 nm 	 Secondary Y Trefoil (5, -3)
 Z20:                 31.23 nm 	 X Pentafoil (5, 5)
 Z21:                -24.02 nm 	 Y Pentafoil (5, -5)
 Z22:                 13.79 nm 	 Secondary Spherical (6, 0)

Total RMS: 	 398.5 nm



In [60]:

    
zv.fringe_bar_chart().show()
plt.show()



In [61]:

    
plots.plotImage(algo.Wconverge, "Final wavefront", show=True)



In [62]:

    
plt.imshow(algo.image)
plt.show()



In [ ]:

    
plt.close('all')

Scaling seems wrong so try to calculate offset for MMT same way...



In [31]:

    
pupsize = 135. * 26.
diameter = 6.5 * u.m
radius = diameter / 2.
focal_length = diameter * 5.3
mmt_ang = np.arctan2(radius, focal_length)
offset = 0.5 * pupsize / np.tan(mmt_ang) * u.um
offset









    Out[31]:





$18603 \; \mathrm{\mu m}$

OK, this matches what we're using there.

Use the Z04 term from the CWFS fit and see if we can calculate a focus correction



In [38]:

    
# this calculates the difference in size between in-focus pupil and observed in um
foc_err = (zv['Z04'] / (nmperasec * u.nm / u.arcsec)) * (100 * u.um / u.arcsec)

# use the tangent of the beam angle to convert pupil size difference to shift in focus position at focal plane.
# convert to focus counts and negate to create correction to send to telescope.
foc_corr = -counts_per_mm * foc_err.to(u.mm).value / np.tan(ang)
foc_corr









    Out[38]:





$9.8131033 \; \mathrm{}$



In [64]:

    
from astropy.nddata import CCDData
from artnfocus.imutils import ARTNreduce, sub_background, find_donuts



In [72]:

    
im = ARTNreduce("rts2focus0012.fits")
#im.data = sub_background(im)
norm = wfs_norm(im)
plt.subplot(projection=im.wcs)
plt.imshow(im, origin='lower', norm=norm)
plt.xlabel('RA')
plt.ylabel('Dec')
plt.show()